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Abstract 

We compare the performance of the Kramers Equation Monte Carlo (KMC) 
Algorithm with that of the Hybrid Monte Carlo (HMC) algorithm for numer- 
ical simulations with dynamical Kogut-Susskind fermions. Using the lattice 
Gross-Neveu model in 2 space-time dimensions, we calculate the integrated 
autocorrelation time of different observables at a number of couplings in the 
scaling region on 16 2 and 32 2 lattices while varying the parameters of the 
algorithms for optimal performance. In our investigation the performance of 
KMC is always significantly below than that of HMC for the observables used. 
We also stress the importance of having a large number of configurations for 
the accurate estimation of the integrated autocorrelation time. 
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Introduction. Numerical simulation is employed as a significant tool in the nonperturba- 
tive investigation of lattice gauge theories. Most of the serious numerical simulations with 
dynamical lattice fermions have so far been carried out with the so-called Hybrid Monte 
Carlo (HMC) algorithm [0. It has several virtues: it is an exact algorithm, easily imple- 
mentable and reasonably efficient with a dynamical critical exponent z — 1 in the large 
trajectory limit. However, with growing demand of more realistic computations, especially 
in QCD, there is a genuine need to search for better algorithms. There is also some evidence 
PI! of lack of reversibility in the discretized classical equations of motion in the molecular 
dynamics part of the HMC. This will make the algorithm inexact. Since this will presumably 
be worse for larger lattice volumes, it is an unwelcome situation for future computations. 

The main alternatives to HMC, available at the moment, are the local bosonic algorithm 
by Luescher [|J and the Kramers Equation Monte Carlo (KMC) MM, first introduced by 
Horowitz ||. In this letter, we concentrate on the KMC algorithm which is a generalized 
version of the HMC algorithm and is exact. The guiding trajectory of KMC includes a 
dissipative term. The optimally tuned KMC also maintains a dynamical critical exponent 
z — 1, but unlike HMC, for arbitrarily short trajectory length, at least for free field theories. 
In practice, the trajectory length is taken to be a single time-step and this presumably would 
help suppress the effects of the possible lack of reversibility. 

For this reason, the KMC certainly merits serious investigation and has received some 
attention recently p],|5],[7j where KMC has been applied to QCD with the gauge group SU(2) 
and to 1+1 dimensional Gross-Neveu model. The general result of these investigations is 
that the KMC has performed comparably with the HMC, although neither algorithm could 
be called clearly more efficient. Apart from the preliminary report M, the studies used 
Wilson lattice fermions. We feel that there is a necessity for more investigation of KMC, 
especially with Kogut-Susskind fermions. 

In the following we describe our investigation of the KMC algorithm using Kogut-Susskind 
or staggered lattice fermions on the 1+1 dimensional Gross-Neveu model. 

The Gross-Neveu model. Since simulating QCD is computationally demanding, we have 
tested the algorithms on a simple non-gauge model which shares some properties of QCD. 
The Gross-Neveu (GN) model [|[ in 2 space-time dimensions is asymptotically free and 
displays chiral symmetry breaking. Hopefully, our experience with the GN model will serve 
us in good stead when confronted with the real problem, QCD. 

The action of the GN model on a 2-dimensional euclidean lattice can be written down 

as: 

s = 7^ E °\ x ) + E Xf(x)A xyXf (y) (1) 

with, 

A xy = o-{x)5 xy + \ E - <W-a»}] (2) 

= 1, Ux) = (3) 

where x an d X are the single- component Kogut-Susskind fermion fields; o is the scalar 
auxiliary field; x and y are lattice sites and / is a flavor index. If rif represents the number 
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of fermion flavors on the lattice, the number of flavors in the continuum limit is 2n/. The 
lattice spacing is taken to be unity. We note that the fermion matrix A xy is real and 
independent of flavor. 

The above action is the most naive transcription of the GN model on a lattice. No 
improvement to reduce lattice artifacts has been used and the a-term in the fermion matrix 
A xy is taken to be single-site. However, this should be enough for our purpose. 



Implementation of the algorithms. We start by writing the classical Hamiltonian 



sites 



1 - 2 + V+5>/(^r 1$ / 



2g 2 2 



(4) 



which evolves the system in the fictitious time r according to the Hamilton's equations in 
the HMC algorithm. The $ fields above are the usual pseudofermionic fields, p's are the 
momenta conjugate to a and are initially obtained from respective gaussian distributions as 
suggested by the quadratic terms in $ and p in H. The initial configuration of a and p are 
evolved in r according to the Hamilton's equations of motion in HMC while $'s are kept 
fixed. The Hamilton's equations are integrated in discrete steps of 5t using the leapfrog 
scheme for N m d steps to complete a molecular dynamics trajectory at the end of which an 
accept/reject step is performed. This is well known and the details of the HMC algorithm 
can be found in many excellent reviews and articles. 

The KMC involves the modification where a term proportional to momenta p is added 
to the momenta itself in each trajectory, 



7T< = e 



+ VI - e-^i (5) 



where rji is a Gaussian noise with zero mean and unit variance and 7 is a tunable parameter. 
The subscript i refers to the number of trajectory. The leap-frog integration is performed 
over a single molecular dynamics step (N md = 1): 

7r i ((Jr/2)=7r i (0) + yF(«7 i (0)) J 

<7 i+1 (5r) = o-j(O) + STTTiiSr), 

7r l+1 (5r) = tt^t/2) + ^-F (a i+1 (8r)) , 

with the force at the site x 

fm) = + E K(^/W + e}(*)n/0»0] » ( 6 ) 

9 f 

where Q and £ are given by, (A^A)£f = flf = A£f. 

After this a Metropolis accept/reject test with acceptance probability P(ai,7Ti — > 
Ci+ijTTi+i) = ™{l,e fl ' l7 "' ri '" f ''''' +1 ' I ' +1 '} is performed. If the new configuration (cr i+1 , 7r i+1 ) 
is accepted, we start again from Eq. || with p i+i = 7Tj + i, and if rejected, we set a i+ i = 
Pi-, Pi+i — ~ 7[ i an d try again starting from Eq. ^. It is important to negate the mo- 
menta in case of rejection to guarantee exactness of the algorithm. Momenta p and the 
pseudofermionic fields $ are refreshed after every fc-trajectories. 



3 



Here 7 and k are tunable parameters. In the limit 7 = 00 and k = 1 KMC reduces to 
HMC algorithm with N md = 1. 

Autocorrelation time. A good measure of how correlated are the configurations generated 
by HMC and KMC, for a particular observable, is the integrated autocorrelation time for 
that observable. For either algorithm, we have measured it for o and the fermion propagator. 

Using the window method suggested by Sokal ||, an estimator for the integrated auto- 
correlation time r int of some observable 9 can be expressed as 

T int(T cu t) — — + E Cg(0) ' ^ 

where Cg(t), the unnormalized autocorrelation function, can be estimated by the following 
expression: 

1 

Ce(t) = —-Y,&-°m + t-9). (8) 

n \ l \ i=i 

9i is the average value of 9 for the i-th configuration and 9 is obtained by averaging over the 
9iS. n is the total number of configurations, a large number, but obviously finite. 

In the above, the factor of 1/2 is purely a convention and T cut is a cut-off for the sum, to 
be chosen judiciously so that one strikes a balance between noise (in case T cut is too large) 
and bias (in case T cut is too small) in the estimation of Tint. in short, the recommendation 
of ref. |9j is to take T cut > CT int , where c is a number between 4 and 10 depending on the 
system and n ^> r int . 

If one plots r int as a function of T cut for a given observable, one can read off the value of 
Ti n t from a plateau region where Tj nt is independent of T cut . Such plateau regions are indeed 
obtained and will be discussed when results are presented in the following. 

Performance tests. The numerical simulations were carried out on 16 2 and 32 2 lattices 
using rif = 4. For HMC, we maintained N mi ib~T ~ 1 and the ansatz for the input vector 
for the matrix inversion at each time-step within the molecular dynamics trajectory was 
simply taken to be the result-vector of the matrix-inversion of the previous time-step. For 
both algorithms, the parameters of the guiding Hamiltonian were not tuned for optimal 
performance and the conjugate gradient inverter was used for fermion matrix inversions 
without any preconditioning. Preconditioning does not seem to work anyway for Kogut- 
Susskind fermions. The parameters of the algorithms were always tuned so as to keep the 
acceptance rate above 90%. 

We tried both cold and hot starts, but for the production runs we always used cold starts 
for both algorithms. At each set of parameters we typically used 1000 iterations in HMC and 
16000 iterations in KMC for equilibration. As we will see, these numbers are substantially 
larger than 20r in4 , suggested by 0, for the observables of interest. 

We determined the scaling region by plotting (a) as a function of 1/g 2 . The scaling 
window is between g = 0.55 — 0.45 for 16 2 lattice and gets slightly beyond 0.45 for 32 2 . For 
the investigation of autocorrelation times, we collected data at various values of the coupling 
g always remaining within the scaling window. 
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For HMC, we varied St from 0.033 to 0.05 always maintaining N m dST ~ 1. For each 
coupling we tried to find the best value for the ratio r = a/2Tj nt where a is the acceptance 
rate. The bigger the ratio, the better is the performance. On 16 2 lattice at g — 0.5 this is 
obtained at St = 0.04 for the observable (a) and roughly also for the fundamental fermion 
propagator. In Fig. la we have plotted T int versus T cut for each of these observables for 16 2 
lattice. The lower curve is for the fermion propagator, and to be specific, we have presented 
the autocorrelation data for zero-spatial-momentum propagator at time separation of 3. 
For the subsequent plots of fermion propagator autocorrelation, we will use this particular 
propagator. For a given T cut , we have measured T int from each of 15 segments, each segment- 
size being 1500 configurations for (a) and 1000 for the fermion propagator. The errors on 
the Tmt data represents the statistical error from the 15 values for each T cut . In subsequent 
plots the same procedure will be followed, although the segment-size and the number of 
segments will be changed appropriately. From Fig. la, one can easily find a plateau region 
from which T int can be read off. One sees also that the fermion propagator has significantly 
smaller integrated autocorrelation time. 

For a reasonably accurate determination of Ti nt we found it absolutely essential to have 
large segment sizes (~ 1000T int ). With smaller segment sizes we have found it difficult to 
find the plateau because of excessive noise. As a general observation, we have noticed that 
the value of T int , if at all one can read it off, is always lower when the segment-size is small. 
In HMC, for both 16 2 and 32 2 lattices, we have actually determined T int at least for (a) for 
a host of other values of g — 0.52, 0.48, 0.46 with segment-size of 100 configurations and 
have rough estimation of autocorrelations at those couplings too. For KMC, the ill-effects 
of having small segment sizes were more dramatic. 

For KMC, we varied St from 0.033 to 0.05 on 16 2 and 32 2 lattices and also varied the 
parameters 7 and k. Within our range of the parameters the near-optimal choice for the 
ratio r at g — 0.5 was St = 0.04, 7 = 1.5, k = 4 on 16 2 lattice and in Fig. lb we present the 
corresponding autocorrelation time data (later, we found that by increasing k to 8, it was 
possible to enhance the performance of KMC somewhat, and thus our choice of parameters 
is not quite optimal). Large statistics was required to produce the figure. The segment-size 
was 72,000 configurations for fermion propagators and 100,000 for (<r), with 5 such segments 
for each observable. Once again, nice plateau regions were found for T int . The propagators 
are again less correlated, but the values of Ti nt for either observable in the plateau region 
are orders of magnitude higher than those obtained in Fig. la for HMC. The acceptance 
rates are usually always higher for KMC, but the huge discrepancy in the autocorrelations, 
even for observables like the propagators, makes the HMC significantly better in effective 
performance, in spite of HMC's higher overhead per trajectory (see Table 1 below). 

We made a qualitative estimate of the higher overhead of HMC over KMC per trajectory 
for all our runs. Ratio of total number of conjugate gradient iterations for one full trajectory 
of HMC to that for one trajectory of KMC (including Hamiltonian calculations in each case) 
gives an accurate estimate of the overhead. For the runs in Figs. la and b, this overhead is 
14.0 with 10% uncertainty. For all our runs at different couplings, lattice sizes and algorithm 
parameters, the HMC overhead was between 10 and 15. 

To illustrate the point of having large enough segment sizes, we plot in Fig. 2 again Ti nt 
versus T cut for (a) on 16 2 lattice at g — 0.5 using the same set of parameters, but this time 
using 6 segments each of size 1000 configurations only. Please note the r inr scale of Fig. 2 as 
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compared to that of Fig. lb. There is only a rough hint of a plateau in Fig.2 and the data is 
very noisy. Also, the actual value of r, n t at the anticipated plateau region of Fig.2 is about 
half of that in Fig. lb. 

In KMC we have investigated integrated autocorrelation times at various other couplings, 
viz., g = 0.52, 0.48, 0.46 with 7 = 1.5, 5r = 0.04 and k = 4 on both 16 2 and 32 2 lattices, 
however, with smaller statistics as in Fig.2 and could estimate Tint for the observables only 
very roughly. Still in each case the autocorrelation time was found to be orders of magnitude 
higher than the corresponding values with HMC, as shown in Table 2 presented later. 

On 16 2 lattice, with k = 4 and 5r = 0.04 we have also varied the parameter 7 from 0.5 to 
2.0 in intervals of 0.5 for a number couplings g from 0.52 to 0.46. The autocorrelation time 
determinations could only be done with relatively small statistics. The optimal 7 was found 
in each case to be qualitatively consistent with the free theory estimate of l/m s where m s 
is the smallest mass of the theory. As a result the optimal 7 shifted towards higher values 
as the coupling g was lowered with all other parameters kept fixed. 

In addition to k = 4, we have performed similar autocorrelation investigations at k = 8 
on 16 2 lattice using 7 = 1.5, 8t — 0.04 at couplings g = 0.50, 0.46. By increasing the 
value of k, there was an indication in our data of an enhancement of the ratio r at the lower 
coupling. Still it was too little to boost it up on par with HMC. 

In Table 1 we present the values of r int and the ratio r for (a) and the fermion propagator 
at g = 0.5 on 16 2 lattice in both HMC and KMC for the data shown in Fig. la and Fig.lb. 
Put together with higher overhead for HMC (dicussed earlier), table 1 then indicates that 
HMC is about 9.5 and 7.5 times more efficient than KMC respectively for ((f)) and for fermion 
propagator. Table 2 collects some results from small statistics runs on both 16 2 and 32 2 
lattices at a few other couplings in the scaling region. The data at g = 0.46 on 16 2 lattice 
has higher statistics than the rest in Table 2, as evident from smaller errors. The KMC data 
on 32 2 lattice are the most inaccurate and the values of Ti nt are likely to be significantly 
below their actual values. 

Conclusions. In our investigation of relative performances of HMC and KMC algorithms 
with dynamical Kogut-Susskind fermions on a simple 2-dimensional Gross- Neveu model, we 
found that the HMC was significantly more efficient in producing uncorrelated configurations 
than the KMC. When other similar investigations 0|| with Wilson fermions did not find a 
significant difference of performance between the two algorithms, our results indicate, taking 
into account the higher overhead per trajectory of the HMC which in our case was always 
between 10 and 15, a difference of about an order of magnitude for all combinations of 
parameters and couplings we tried. How much of our conclusion would carry over to QCD 
with Kogut-Susskind fermions is an open question. 

Our choice of parameters for KMC might not have been optimal in all cases, but it 
becomes apparent that KMC requires a lot of tuning and our investigation shows that we 
were unable to find a set of parameters which could make KMC even remotely competitive. 

We emphasize again the importance of having a large number of configurations for the 
estimation of T int . Without it, the data is not just noisy, the roughly estimated T int is well 
below the actual value. 
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FIGURE CAPTIONS 

Fig. 1. Integrated autocorrelation times for (a) and fermion propagator (zero spatial 
momentum with time-separation 3) at g — 0.5 on 16 2 lattice for (a) HMC with 5t = 0.04, 
(b) KMC with 5t = 0.04, 7 = 1.5, k = 4. 



Fig. 2. Low statistics data on r int for (a) in KMC at g = 0.5 on 16 2 lattice with the same 
set of parameters as in Fig. lb. 

TABLE CAPTIONS 

Table 1. r int and r for g = 0.50 obtained on 16 2 lattice (HMC: 5r = 0.04; KMC: 5r = 0.04, 
7 = 1.5, k = 4). 



Table 2. r int and r for various g obtained on 16 2 and 32 2 lattice for (a) (HMC: 5r = 0.04; 
KMC: 5t = 0.04, 7 = 1.5, k = 4). 
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TABLES 



Algorithm 




propagator 




Tint 


r 


Tint 


r 


HMC 


1.72(28) 


29(5) 


0.83(11) 


63(8) 


KMC 


228(27) 


0.22(2) 


83(6) 


0.60(5) 



TABLE 1. 





g = 0.52 


g = 0.48 


g = 0.46 


Tint 


r 


Tint 


r 


Tint 


r 


16 2 


HMC 
KMC 


1.22(48) 
128(64) 


45(18) 
0.50(25) 


1.30(45) 
108(26) 


40(13) 
0.47(11) 


1.05(20) 
138(23) 


48(9) 
0.36(6) 


32 2 


HMC 
KMC 


1.47(38) 
124(97) 


33(8) 
0.84(66) 


1.47(50) 
142(107) 


33(11) 
0.63(47) 


1.56(31) 


30(6) 



TABLE 2. 
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FIGURES 




FIG. 2. 



